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Abstract 

In developing mathematical models of systems from a given input“Output 
data sequence, the choice of the sampling interval and the selection of the 
order of the model in time-series analysis pose difficult problems. Band- 
limited (up to 15 Hz) random torque perturbations were applied to the human 
ankle joint. The applied torque input, the angular rotation output, and the 
electromyographic activity using surface electrodes from the extensor and flexor 
muscles of the ankle joint were recorded. Autoregressive moving average models 
were developed. A parameter constraining technique is applied to develop more 
reliable models. It is shown that the asymptotic behavior of the system must 
be taken into account during parameter optimization to develop predictive models. 

INTRODUCTION 

In a series of previous papers (Agarwal and Gottlieb, 1977 a,b; Gottlieb 
and Agarwal, 1978; Gottlieb, Agarwal, and Penn, 1978) we have attemped to 
describe quantitatively the neuromuscluar system dynamics to applied sinusoidal 
and band-limited gaussian torque perturbations. In these studies, the compliance 
of the joint was calculated using Fourier series analysis for sinusoidal and 
power spectral density methods for random perturbations. Although linear 
analysis methods were used, the system is known to be nonlinear and the parameter 
values such as the joint viscous and stiffnegs coefficients are functions of 
the level of neuromusclar activity. 

The purpose of the present paper is to apply time series analysis methods 
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to study the input-output behavior of the neuromuscular system. The time series 
method is very parsimonious in the use of parameters to represent the model 
structure. Normalized residual criterion (NRC) will be used to estimate the 
model order (For details of this method see Suen and Liu, 1977; Osaf o-Charles 
et. al., 1980), 

Our previous analysis was limited to analysis of the angular rotation data 
and calculation of joint compliance. The electromyographic (EMG) data was not 
analysed due to inherent difficulties in representing this output by linear 
transfer functions. The time series approach allows nonlinear representations 
as long as the model is linear in parameter space. 

Dufresne, Soechting and Terzuolo (1978) used pseudo-random torque pulses 
to study the human forearm response. They developed a model of the EMG in 

terms of the limb position and its derivatives in the following form: 

• ( 1 ) 
EMG(t) = A e(t - d) + B 0(t - d) -H C 0(t - d) ^ ' 

where A, B, and C are constant parameters and d is the time delay. They 
found that the motor output depends primarily on the angular velocity of the 
joint. The time delay was found to be about A7 msec. 

In a subsquent study, Dufresne, Soechting, and Terzuolo (1979) used 
different time delay parameters for position and its two derivatives. The 
best estimates for the time delays were found to be 86 msec for position, 25 msec 
for velocity, and A5 msec for acceleration. The physiological processes 
associated with these varying delays are not clear. Soechting and Dufresne 
(1980) found that the linear model given in equation (1) predicted 80% of the 
EMG response. 

Our analysis of the EMG using time series shows that the autoregressive 
terms of the EMG are important and cannot be Ignored as was done in the Dufresne 
et. al (1978) model. 
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METHODS 


These experiments were done using normal human subjects. A subject sat 
in a chair with the right foot strapped to a footplate which could rotate about 
a horizontal, dorsal-plantar axis through the medial malleolus. The plate could 
be rotated by a DC torque motor. A band-limited gaussian ( 0-15 Hz ) signal was 
prerecorded from a noise generator. These time-varying signals were superimposed 
on a biasing mean motor torque level. The subject was instructed to try to 
maintain a constant mean force against the bias torque of the motor so that the 
ankle joint movement was nearly symmetrical with respect to the reference angle- 
The input was applied for 30 sec or more and the data continuously recorded on 
a digital tape. 

The torque was measured by a strain gauge bridge on the side arms of the 
footplate. Angular rotation was measured by a continuous capacitive transducer. 
The EMGs were recorded from disc surface electrodes taped over the bellies of 
the soleus (SM) and the anterior tibial (TA) muscles. These were amplified 
full-wave rectified and passed through an averaging filter (10 msec averaging 
time) bofore recording. A computer generated the motor drive voltage at 
a conversion rate of 250/s and digitized data on four input channels. The angle 
and the torque signals were sampled at a rate of 250/s and filtered EMGs at 
a rate of 500/s. The data analysis was done off-line using the Minitab 2 
statistical software package on an IBM 370 computer. 
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The Normalized Residual Criterion 


Time-series analysis can be extended to obtain discrete linear transfer 
functions of systems having an input x(t) and output y(t). By x(t) and y(t) we 
mean pairs of observations that are available at equispaced Intervals of time. 

The behavior of the dynamic system can be adequately represented by the present 
and past responses and the current and past Inputs of the systems. We denote this 
process as transfer function (TF) models (n,m) and write its equation as 


y(t) = ao + '*iy(t - 1) + Oj^yCt - n) + Sox(t)+. . .+6^^ x(t-n)+v(t) 

In (2) the parameters to be estimated are ao>...«a , « n, and 

n m 

ra. The time seiies v(t) is a random term measuring the difference between the 
response y(t) and the variables used to expalin the time-series data. The 
parameter ao measures the mean output. 

Equation (2) reduces to an autoregressive model (AR(n)) if x(t) is omitted 
from the model, and reduces to a moving average model (MA(m)) if lags of y are 
omitted. The following assumptions will be made concerning v(t) for a given 
ouput time sequence y(t), t = [0,T], 


1) 

E[v(t)] = 0 


:) 

E['i(i)v( j) ] = 


where 
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Define t - 1, 2, •••,T-n. (4) 
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Note that in the discussion below V and Y are vectors such that 



■ v(l) 


y(l) ^ 



v(2) 

• 

• 

Y* 

y(2) 



• 

_v(T-n) J 


..y(T-n) _ 

(6) 


(7) 


Squaring (4) and normalizing by the total sum of squares, we have 

JMli . . | |y-»n-ra,y - r3j» || ^ . __ ,, 

and therefore 

e[ ||v|| 2 ]= E[|lYil2e(n,m.T)] . (8) 

Since y(t), the data series, is deterministic, (8) can be rewritten as 

[ llvIP ]- IIyM^E [e(n, m, T) ]. 


E|_ llVlP 

From (6) we have 


e[ I lv| |2 ]= E [v2(D +v2(2) + **•+ v2(T-n)] , 
\;hich by assumption 2) in (3) reduces to 

E [imp] ■ [y-njoj . 

Substitution of (11) in (9) we have 


(9) 


( 10 ) 


( 11 ) 


l^c (n,m,T) j 




( 12 ) 


and bs’ assumption 3), (12) becomes 


|^c(n, m, T)j= 


To! 


(13) 


Y 


The quantity c(n,ra,T) depends cn n, m, and T and is proportional to the 
nortmlized variance of the regression for a given n and m. If this ratio is 
mininized over n and m, then the data fit as measured by the correlation coeffi- 
cient P will be maximized. Note that 
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p 


1 - 


2 


1/2 



( 14 ) 


( 15 ) 


where T, being a constant for the data, Is omitted in the optimization procedure, 
and £(n,m) is the minimum value for e(n,m). This optimization technique is the 
so called the Normalized Residual Criterion. 
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RESULTS 


The mathematical modeling problem was considered in two separate parts. 

For the first model the applied torque is the input and the resulting angular 
rotation of the joint is the output of the system. For the second model, the 
angular rotation (and its derivatives) is considered as the input to the 
system and the resulting stretch reflex electromyographic activity is consi- 
dered as the output. It should be emphasized that the angular rotation is 
the net result of two torque inputs applied at the joint; one by the external 
motor torque and the other muscle forces produced by the stretch reflex mech- 
anism. These mechanisms are also responsible for a significant contribution 
to the joint viscous and elastic properties. Figure 1 shows a sample of the 
data at 4 msec sampling interval. The velocity was obtained by digital 
differentiation. 

Angle-Torque Model 

Although the data was recorded for 30 seconds at each input (Agarwal and 
Gottlieb, 1977b), this method does not require such long data records which 
would also use too much computer time. The time series analysis was done using 
only two-seconds of the data record. (The first two-seconds of the data were not 
used to allow the turn-on transients to die out). 

The values of f(n,m'i were computed for a given data record and then plotted 
against different values of n (see Figure 2). Tne data sampling interval in 
this case is 4 msec. This analysis clearly indicates that n * 2 and m • 0 is 
adequate to model this data. The same data was analyzed again using the sampling 
intervals of 12, 20, 40, and 60 msec. Figure 3 shows the e(n,m) values for the 
sampling interval of 20 msec. Note that the minimum value of the normalized 
residual is about 60 times of that in the first case. For 40 and 60 msec sampling 
e(n,m) did not reach asNTnptotic values even for model order of (8,8). The norm- 
alized residuali' at 12 and 20 msec sampling implied a model order of (3,1). 
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Figure 4 shows the actual angular rotation data (2 to 4 sec Interval used 
in this analysis ), the regression fit and the predicted output using 4 msec 
sampling and model order of (3,1). The regression fit is obtained by using the 
equation 

9(t) - ao + ai 0(t - 1) + 02 6 (t - 2) + 03 6 (t - 3) + BoT(t) + BiT(t - 

(16) 

The error between the actual data and the regression fit is nearly i:ero. 

The correlation coefficient is p ■ 0.999. However, when this mode* is sed to 
predict the output using the first three data output values as the initial con- 
ditions, the predicted output is a poor approximation of the actual data (see 
Figure 4). Figure 5 shows the observed angle and the predicted model values for 
model ordtrs of (3,1), (7,1), (9,1), and (14,1). Even the fourteenth order 
mv>del is not able to adequately reproduce the data sequence. These models are 
not able to capture the steady state (or long term) behavior of the system. 
Osafo-Charles, et al., (1980) showed tiiat to develop better predictive models, 
the TF(n,m) models must be constrained to incorporate the steady state response 
of the system. 

C onstrained Model 

Consider the estimated model given by equation (16). Under conditions 
of equilibrium 

0(t) • 6 (t - 1) • 0 (t - 2) « 0(t - 3) - 0 ^ 

e 

and 

T(t) • T<t - 1) - Tp 

wliere 6 ^ and are the steady state response and input respectively. At physical 
and statistical equilibrium, with 6 (t) ■ 6 ^ and T(t) ■ T*. equaition (16) became 


*3^ “ Uj 0 ^, ^2 ®e ^^3 ®e + dj T^, 


(17) 


or 


bo 6i 


1 - a\- c»2 - aj 


K 


(18) 





f 


where (18) expresties the steady state gain In terms of the parameters of the model. 

The value of g was approximated by the slope of the curve of torque vs. angular 
rotation in the relaxed ankle during sinusoidal oscillation at 0.1 (Gottlieb 
and Agarwal, 1978). 

For 6 g ■ g Tg to be true, we must have 


Bq ” g(l - Oj - 02 - O3)- 


From equrtions (19) and (16), we get 

[e ;t) - g T (t)] - oi 

fi(t - 1) - g T(t) 

+ CH2 

0(t - 2) - g T(t) 

+ 03 

e(t - 3) - g T(t) 


T(t - 1) - T (t)] 


(19) 


( 20 ) 


Regression analysis is used again to estimate the parameters a^, 02 , 

03 and for a given value of gain g. Bq is then obtained using (19). Figure 
6 shows the output angle and predicted model response for a constrained model 
with gains of g » 6.5, 7.5, and 8.5. The gain value of 8.5 was considered to 
provide the best fit in terms of the minimum estimated standard deviation of the 
regression. 

The transfer function for the unconstrained model is: 

H(z) - 0.00239 - 0.00024 

-1 -2 ( 21 ) 

1 - 2.678 z + 2.399 2 - 0.7191 z 


For the constrained model with a slope of 8.5, the transfer function is: 

H(z) - 0.00238 + 0.00017 ^^ 2 ) 

1 - 2.731 z"* + 2.503 z“" - 0.7717 z'^ 


L 
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EMC Model 


Our efforts to taodel EMC as a function of either the angular rotation 
or the velocity or a combination of both were not successful. As was noted 
by Dufresne et al. (1978), the velocity of rotation is the most significant 
input due to spindle properties (Matthews, 1972). However only those 
components of velocity which stretch the spindle contribute to the EMG of 
the stretched muscle. (The spindle is silent during shortening). Therefore, 
a new velocity signal represneting only the stretching velocities was 
dt. fined as: 


• « * 

(t) - 0 (t) if 6^0 

- 0 if 6 < 0 


(23) 


The normalized residual analysis indicated a model order of (A,l) using 
soleus LMG as the ouput and 6^ as the input signal. The predicted output of 
the unconstrained model and its comparison with the actual EMG signal is 
shown in Figure 7. Since the EMG signal is a full-wave rectified and filtered 
(using an averaging filter) signal, it has only negative values (because of 
negative filter gain). The predicted value of EMG is a poor approxim ition of 
the data. 

A const! aiued model was developed ui»ing a similar approach as outlined 
earlier. Figure (8) shows tie predicted EMC and the actual data at three 
values the gein parameter. The gain of -0.003 was considered to give the 

nost appropriate fit. 

For the constrained model with a slope of -0.005, the transfer function 


is; 


H(z) 


EMG 


-.004164 .00314 z~ ^ ( 24 ) 

1 - 1.19 z~' + .6685 z"- - .2947 z"^ + .02104 z’^* 
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CONCLUSIONS 


The time series approach is a powerful and versatile technique in 
developing time domain models from a given input-output data sequence. Norm- 
alized residual criterion allows effective prediction of the model order. 
Models developed in this manner may be satisfactory, but may not be good 
predicative models. It is recommended that constrained parameter modeling 
which allows incorporating the steady-state behavior be used to obtain 
better predictive models. 
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